dat <- readxl::read_xlsx("/Volumes/manhnd/Yi_roe_data.xlsx")
# filter variables for analyis
dat <- dat %>% dplyr::select("serialno", "country", "papm", "gender", "age", "education", "marital", "occupation", "live_area", "mobile", "internet")
# Remove the gender as others (18/2912)
dat1 <- dat %>%
filter(gender != 3 & gender != 4) %>%
mutate(marital = dplyr::recode(marital, `4` = 3, `5` = 3))
# Transform the valuables
dat1 <- dat1 %>%
mutate(
papm = factor(papm, ordered = TRUE),
gender = factor(gender, labels = c("Male", "Female")),
marital = factor(marital, labels = c("Never married", "Married/living with a partner", "Separated/widowed/divorced")),
education = factor(education, labels = c("No formal education", "Primary school", "Secondary school", "Diploma", "Degree")),
occupation = factor(occupation, labels = c("Employed", "Unemployed", "Student", "Retiree", "Homemaker")),
live_area = factor(live_area, labels = c("Urban", "Peri-urban", "Rural", "Slums")),
mobile = factor(mobile, labels = c("No", "Yes feature phone", "Yes smart phone")),
internet = factor(internet, labels = c("No", "Yes")),
country = factor(country)
)
summary(dat1)
## serialno country papm gender age
## Min. : 1.0 Philippines:459 1:1596 Male :1365 Min. :18.00
## 1st Qu.: 726.2 Uganda :395 2: 291 Female:1529 1st Qu.:22.00
## Median :1453.5 Indonesia :368 3: 285 Median :30.00
## Mean :1455.6 Zimbabwe :292 4: 179 Mean :34.95
## 3rd Qu.:2185.8 Cameroon :289 5: 543 3rd Qu.:44.00
## Max. :2912.0 Bangladesh :280 Max. :90.00
## (Other) :811
## education marital
## No formal education: 173 Never married :1221
## Primary school : 476 Married/living with a partner:1395
## Secondary school :1182 Separated/widowed/divorced : 278
## Diploma : 405
## Degree : 658
##
##
## occupation live_area mobile internet
## Employed :1294 Urban :1066 No : 257 No : 839
## Unemployed: 566 Peri-urban: 374 Yes feature phone: 626 Yes:2055
## Student : 681 Rural :1281 Yes smart phone :2011
## Retiree : 65 Slums : 173
## Homemaker : 288
##
##
str(dat1)
## tibble [2,894 × 11] (S3: tbl_df/tbl/data.frame)
## $ serialno : num [1:2894] 1 2 3 4 5 6 7 8 9 10 ...
## $ country : Factor w/ 9 levels "Bangladesh","Cameroon",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ papm : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 1 1 1 1 1 1 1 1 1 1 ...
## $ gender : Factor w/ 2 levels "Male","Female": 2 2 2 2 2 1 1 2 1 2 ...
## $ age : num [1:2894] 29 36 52 27 31 31 70 50 49 45 ...
## $ education : Factor w/ 5 levels "No formal education",..: 2 1 1 2 2 2 1 1 1 1 ...
## $ marital : Factor w/ 3 levels "Never married",..: 2 2 3 2 2 2 2 2 2 2 ...
## $ occupation: Factor w/ 5 levels "Employed","Unemployed",..: 5 5 1 5 5 1 2 5 1 5 ...
## $ live_area : Factor w/ 4 levels "Urban","Peri-urban",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ mobile : Factor w/ 3 levels "No","Yes feature phone",..: 2 2 2 2 2 2 2 1 2 2 ...
## $ internet : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
dat2 <- dat1 %>%
select("country", "papm") %>%
group_by(country, papm) %>%
summarise(n = n()) %>%
mutate(proportion_papm = n / sum(n)) %>%
ungroup()
## `summarise()` has grouped output by 'country'. You can override using the
## `.groups` argument.
ggplot(data = dat2, aes(x = country, y = proportion_papm)) +
geom_col(width = 0.5) +
facet_grid(~papm) +
coord_flip()
dat2 <- dat %>%
select("country", "education") %>%
group_by(country, education) %>%
summarise(n = n()) %>%
mutate(proportion_education = n / sum(n)) %>%
ungroup()
## `summarise()` has grouped output by 'country'. You can override using the
## `.groups` argument.
ggplot(data = dat2, aes(x = country, y = proportion_education)) +
geom_col(width = 0.5) +
facet_grid(~education) +
coord_flip()
m1 <- polr(formula = papm ~ country + gender + age + education + marital + occupation + live_area + mobile + internet, data = dat1, Hess = TRUE)
summary(m1)
## Call:
## polr(formula = papm ~ country + gender + age + education + marital +
## occupation + live_area + mobile + internet, data = dat1,
## Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## countryCameroon 1.030800 0.304414 3.38618
## countryIndia 3.224494 0.272151 11.84819
## countryIndonesia 1.723045 0.279535 6.16397
## countryKenya 2.155655 0.310191 6.94944
## countryNepal 3.933541 0.282818 13.90837
## countryPhilippines 1.061534 0.286374 3.70681
## countryUganda 4.131868 0.272311 15.17337
## countryZimbabwe 4.363371 0.287158 15.19501
## genderFemale -0.005772 0.085817 -0.06726
## age -0.001336 0.003866 -0.34557
## educationPrimary school 0.480407 0.190301 2.52446
## educationSecondary school 0.858222 0.195790 4.38338
## educationDiploma 1.083491 0.222924 4.86036
## educationDegree 1.204571 0.220266 5.46872
## maritalMarried/living with a partner 0.243062 0.122772 1.97978
## maritalSeparated/widowed/divorced 0.602561 0.181871 3.31312
## occupationUnemployed -0.183360 0.118656 -1.54530
## occupationStudent -0.156333 0.134595 -1.16151
## occupationRetiree -0.556455 0.284433 -1.95637
## occupationHomemaker -0.404096 0.151572 -2.66603
## live_areaPeri-urban -0.019113 0.135460 -0.14110
## live_areaRural 0.056872 0.116823 0.48682
## live_areaSlums 0.341699 0.225495 1.51533
## mobileYes feature phone 0.576591 0.161966 3.55994
## mobileYes smart phone 0.217389 0.212906 1.02106
## internetYes 0.491155 0.168225 2.91964
##
## Intercepts:
## Value Std. Error t value
## 1|2 4.2778 0.3845 11.1264
## 2|3 4.9192 0.3873 12.7018
## 3|4 5.5729 0.3902 14.2807
## 4|5 6.0354 0.3921 15.3938
##
## Residual Deviance: 6159.317
## AIC: 6219.317
# First fit the model with interaction between country and all of the independent variables but the model cannot run. Then I removed the interaction between live_area and country, model can run as below
m2 <- polr(formula = papm ~ country + gender + age + education + marital + occupation + live_area + mobile + internet + gender:country + age:country + marital:country + occupation:country + mobile:country + internet:country + education:country, data = dat1, Hess = TRUE)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(m2)
## Call:
## polr(formula = papm ~ country + gender + age + education + marital +
## occupation + live_area + mobile + internet + gender:country +
## age:country + marital:country + occupation:country + mobile:country +
## internet:country + education:country, data = dat1, Hess = TRUE)
##
## Coefficients:
## Value Std. Error
## countryCameroon 14.609312 1.400e+00
## countryIndia 24.913158 1.278e+00
## countryIndonesia -3.055593 1.099e+00
## countryKenya 26.698692 1.697e+00
## countryNepal 27.589432 1.359e+00
## countryPhilippines 23.911632 1.497e+00
## countryUganda 26.874774 1.138e+00
## countryZimbabwe 16.276852 1.281e+00
## genderFemale 0.354487 5.578e-01
## age 0.077100 3.809e-02
## educationPrimary school -4.080201 3.335e-01
## educationSecondary school 9.866479 7.889e-01
## educationDiploma 11.217361 6.772e-01
## educationDegree 11.537265 5.694e-01
## maritalMarried/living with a partner -0.857760 1.184e+00
## maritalSeparated/widowed/divorced -2.453037 2.089e-01
## occupationUnemployed -16.542871 1.461e-01
## occupationStudent 0.048204 1.308e+00
## occupationRetiree -18.167315 3.829e-01
## occupationHomemaker -12.483776 2.274e-01
## live_areaPeri-urban -0.090668 1.443e-01
## live_areaRural -0.001955 1.294e-01
## live_areaSlums 0.249086 2.547e-01
## mobileYes feature phone 11.977607 8.830e-01
## mobileYes smart phone 11.867155 1.091e+00
## internetYes 0.097749 1.624e+00
## countryCameroon:genderFemale -0.475020 6.699e-01
## countryIndia:genderFemale 0.074489 6.111e-01
## countryIndonesia:genderFemale -0.611408 6.207e-01
## countryKenya:genderFemale -0.179755 6.525e-01
## countryNepal:genderFemale -0.674883 6.148e-01
## countryPhilippines:genderFemale -0.091360 6.396e-01
## countryUganda:genderFemale -0.477477 5.858e-01
## countryZimbabwe:genderFemale -0.918931 6.375e-01
## countryCameroon:age -0.074785 4.493e-02
## countryIndia:age -0.063382 3.976e-02
## countryIndonesia:age -0.033517 4.247e-02
## countryKenya:age -0.031155 4.203e-02
## countryNepal:age -0.064597 4.016e-02
## countryPhilippines:age -0.080832 4.119e-02
## countryUganda:age -0.073202 3.865e-02
## countryZimbabwe:age -0.122026 4.188e-02
## countryCameroon:maritalMarried/living with a partner 1.420520 1.273e+00
## countryIndia:maritalMarried/living with a partner 0.976123 1.236e+00
## countryIndonesia:maritalMarried/living with a partner 0.718878 1.232e+00
## countryKenya:maritalMarried/living with a partner 0.783927 1.265e+00
## countryNepal:maritalMarried/living with a partner 0.381139 1.243e+00
## countryPhilippines:maritalMarried/living with a partner 1.822441 1.282e+00
## countryUganda:maritalMarried/living with a partner 0.919528 1.221e+00
## countryZimbabwe:maritalMarried/living with a partner 2.381148 1.259e+00
## countryCameroon:maritalSeparated/widowed/divorced -12.787097 8.655e-08
## countryIndia:maritalSeparated/widowed/divorced 2.628711 6.726e-01
## countryIndonesia:maritalSeparated/widowed/divorced 0.975907 6.736e-01
## countryKenya:maritalSeparated/widowed/divorced 2.709081 5.056e-01
## countryNepal:maritalSeparated/widowed/divorced 3.123797 6.262e-01
## countryPhilippines:maritalSeparated/widowed/divorced 3.826500 6.370e-01
## countryUganda:maritalSeparated/widowed/divorced 2.576574 4.045e-01
## countryZimbabwe:maritalSeparated/widowed/divorced 4.439632 5.209e-01
## countryCameroon:occupationUnemployed 16.766464 4.480e-01
## countryIndia:occupationUnemployed 16.584719 4.851e-01
## countryIndonesia:occupationUnemployed 15.559812 4.244e-01
## countryKenya:occupationUnemployed 16.774798 4.632e-01
## countryNepal:occupationUnemployed 14.736578 5.290e-01
## countryPhilippines:occupationUnemployed 16.481324 4.856e-01
## countryUganda:occupationUnemployed 16.816579 2.401e-01
## countryZimbabwe:occupationUnemployed 15.691493 3.303e-01
## countryCameroon:occupationStudent -0.162487 1.393e+00
## countryIndia:occupationStudent 0.259397 1.358e+00
## countryIndonesia:occupationStudent -0.397042 1.361e+00
## countryKenya:occupationStudent -0.375236 1.443e+00
## countryNepal:occupationStudent -0.461873 1.359e+00
## countryPhilippines:occupationStudent 0.709372 1.412e+00
## countryUganda:occupationStudent 0.011918 1.360e+00
## countryZimbabwe:occupationStudent 0.052282 1.394e+00
## countryCameroon:occupationRetiree 0.304758 1.344e-08
## countryIndia:occupationRetiree 17.765439 6.078e-01
## countryIndonesia:occupationRetiree 1.689434 1.391e-07
## countryKenya:occupationRetiree 18.059350 1.465e+00
## countryNepal:occupationRetiree 16.820388 5.936e-01
## countryPhilippines:occupationRetiree 18.007567 8.731e-01
## countryUganda:occupationRetiree 17.472933 7.974e-01
## countryZimbabwe:occupationRetiree 36.326893 6.540e-08
## countryCameroon:occupationHomemaker 11.750979 9.985e-01
## countryIndia:occupationHomemaker 11.743691 3.820e-01
## countryIndonesia:occupationHomemaker 11.075414 6.715e-01
## countryKenya:occupationHomemaker 13.993172 7.158e-01
## countryNepal:occupationHomemaker 11.784830 4.632e-01
## countryPhilippines:occupationHomemaker 12.831394 6.032e-01
## countryUganda:occupationHomemaker 12.899331 3.283e-01
## countryZimbabwe:occupationHomemaker 10.611018 9.371e-01
## countryCameroon:mobileYes feature phone -12.848847 1.267e+00
## countryIndia:mobileYes feature phone -11.411442 9.488e-01
## countryIndonesia:mobileYes feature phone -14.822680 2.206e+00
## countryKenya:mobileYes feature phone -11.090439 1.084e+00
## countryNepal:mobileYes feature phone -12.424780 1.400e+00
## countryPhilippines:mobileYes feature phone -11.984637 1.099e+00
## countryUganda:mobileYes feature phone -11.399280 9.077e-01
## countryZimbabwe:mobileYes feature phone -8.440743 1.168e+00
## countryCameroon:mobileYes smart phone -11.923094 1.571e+00
## countryIndia:mobileYes smart phone -11.281882 1.357e+00
## countryIndonesia:mobileYes smart phone -13.735593 2.200e+00
## countryKenya:mobileYes smart phone -10.639113 1.271e+00
## countryNepal:mobileYes smart phone -12.697334 1.963e+00
## countryPhilippines:mobileYes smart phone -13.246191 1.279e+00
## countryUganda:mobileYes smart phone -11.300390 1.198e+00
## countryZimbabwe:mobileYes smart phone -10.451132 1.268e+00
## countryCameroon:internetYes -0.142093 1.898e+00
## countryIndia:internetYes 0.343617 1.800e+00
## countryIndonesia:internetYes 19.075766 1.099e+00
## countryKenya:internetYes -0.298728 1.683e+00
## countryNepal:internetYes 1.112136 2.262e+00
## countryPhilippines:internetYes 2.030549 1.746e+00
## countryUganda:internetYes 0.590445 1.694e+00
## countryZimbabwe:internetYes -0.705238 1.672e+00
## countryCameroon:educationPrimary school 13.077706 6.222e-01
## countryIndia:educationPrimary school 4.762476 5.766e-01
## countryIndonesia:educationPrimary school 16.204295 7.376e-01
## countryKenya:educationPrimary school 0.185179 1.174e+00
## countryNepal:educationPrimary school 3.920423 5.980e-01
## countryPhilippines:educationPrimary school 3.667566 8.605e-01
## countryUganda:educationPrimary school 4.581924 3.971e-01
## countryZimbabwe:educationPrimary school 17.735081 7.247e-01
## countryCameroon:educationSecondary school 0.958803 8.201e-01
## countryIndia:educationSecondary school -8.921375 9.204e-01
## countryIndonesia:educationSecondary school 1.257021 7.397e-01
## countryKenya:educationSecondary school -12.508199 1.331e+00
## countryNepal:educationSecondary school -9.343525 9.348e-01
## countryPhilippines:educationSecondary school -9.735912 1.081e+00
## countryUganda:educationSecondary school -9.409909 8.247e-01
## countryZimbabwe:educationSecondary school 3.027823 7.846e-01
## countryCameroon:educationDiploma 0.244220 7.472e-01
## countryIndia:educationDiploma -9.866073 8.812e-01
## countryIndonesia:educationDiploma 0.688810 6.993e-01
## countryKenya:educationDiploma -14.272035 1.290e+00
## countryNepal:educationDiploma -9.862474 8.626e-01
## countryPhilippines:educationDiploma -11.230641 1.041e+00
## countryUganda:educationDiploma -10.547329 7.987e-01
## countryZimbabwe:educationDiploma 1.256645 7.247e-01
## countryCameroon:educationDegree -0.181528 6.094e-01
## countryIndia:educationDegree -9.934391 7.724e-01
## countryIndonesia:educationDegree 0.151529 5.295e-01
## countryKenya:educationDegree -13.209072 1.288e+00
## countryNepal:educationDegree -10.477752 7.820e-01
## countryPhilippines:educationDegree -11.690755 9.826e-01
## countryUganda:educationDegree -11.004074 1.124e+00
## countryZimbabwe:educationDegree 0.496835 5.932e-01
## t value
## countryCameroon 1.043e+01
## countryIndia 1.950e+01
## countryIndonesia -2.779e+00
## countryKenya 1.574e+01
## countryNepal 2.030e+01
## countryPhilippines 1.597e+01
## countryUganda 2.362e+01
## countryZimbabwe 1.270e+01
## genderFemale 6.355e-01
## age 2.024e+00
## educationPrimary school -1.224e+01
## educationSecondary school 1.251e+01
## educationDiploma 1.656e+01
## educationDegree 2.026e+01
## maritalMarried/living with a partner -7.245e-01
## maritalSeparated/widowed/divorced -1.174e+01
## occupationUnemployed -1.132e+02
## occupationStudent 3.684e-02
## occupationRetiree -4.744e+01
## occupationHomemaker -5.490e+01
## live_areaPeri-urban -6.285e-01
## live_areaRural -1.511e-02
## live_areaSlums 9.779e-01
## mobileYes feature phone 1.356e+01
## mobileYes smart phone 1.088e+01
## internetYes 6.018e-02
## countryCameroon:genderFemale -7.091e-01
## countryIndia:genderFemale 1.219e-01
## countryIndonesia:genderFemale -9.851e-01
## countryKenya:genderFemale -2.755e-01
## countryNepal:genderFemale -1.098e+00
## countryPhilippines:genderFemale -1.428e-01
## countryUganda:genderFemale -8.151e-01
## countryZimbabwe:genderFemale -1.442e+00
## countryCameroon:age -1.664e+00
## countryIndia:age -1.594e+00
## countryIndonesia:age -7.893e-01
## countryKenya:age -7.413e-01
## countryNepal:age -1.609e+00
## countryPhilippines:age -1.962e+00
## countryUganda:age -1.894e+00
## countryZimbabwe:age -2.914e+00
## countryCameroon:maritalMarried/living with a partner 1.116e+00
## countryIndia:maritalMarried/living with a partner 7.898e-01
## countryIndonesia:maritalMarried/living with a partner 5.834e-01
## countryKenya:maritalMarried/living with a partner 6.198e-01
## countryNepal:maritalMarried/living with a partner 3.066e-01
## countryPhilippines:maritalMarried/living with a partner 1.421e+00
## countryUganda:maritalMarried/living with a partner 7.534e-01
## countryZimbabwe:maritalMarried/living with a partner 1.892e+00
## countryCameroon:maritalSeparated/widowed/divorced -1.477e+08
## countryIndia:maritalSeparated/widowed/divorced 3.908e+00
## countryIndonesia:maritalSeparated/widowed/divorced 1.449e+00
## countryKenya:maritalSeparated/widowed/divorced 5.359e+00
## countryNepal:maritalSeparated/widowed/divorced 4.989e+00
## countryPhilippines:maritalSeparated/widowed/divorced 6.007e+00
## countryUganda:maritalSeparated/widowed/divorced 6.369e+00
## countryZimbabwe:maritalSeparated/widowed/divorced 8.523e+00
## countryCameroon:occupationUnemployed 3.743e+01
## countryIndia:occupationUnemployed 3.419e+01
## countryIndonesia:occupationUnemployed 3.666e+01
## countryKenya:occupationUnemployed 3.622e+01
## countryNepal:occupationUnemployed 2.786e+01
## countryPhilippines:occupationUnemployed 3.394e+01
## countryUganda:occupationUnemployed 7.003e+01
## countryZimbabwe:occupationUnemployed 4.750e+01
## countryCameroon:occupationStudent -1.166e-01
## countryIndia:occupationStudent 1.910e-01
## countryIndonesia:occupationStudent -2.918e-01
## countryKenya:occupationStudent -2.601e-01
## countryNepal:occupationStudent -3.400e-01
## countryPhilippines:occupationStudent 5.024e-01
## countryUganda:occupationStudent 8.762e-03
## countryZimbabwe:occupationStudent 3.751e-02
## countryCameroon:occupationRetiree 2.267e+07
## countryIndia:occupationRetiree 2.923e+01
## countryIndonesia:occupationRetiree 1.215e+07
## countryKenya:occupationRetiree 1.233e+01
## countryNepal:occupationRetiree 2.834e+01
## countryPhilippines:occupationRetiree 2.063e+01
## countryUganda:occupationRetiree 2.191e+01
## countryZimbabwe:occupationRetiree 5.555e+08
## countryCameroon:occupationHomemaker 1.177e+01
## countryIndia:occupationHomemaker 3.075e+01
## countryIndonesia:occupationHomemaker 1.649e+01
## countryKenya:occupationHomemaker 1.955e+01
## countryNepal:occupationHomemaker 2.544e+01
## countryPhilippines:occupationHomemaker 2.127e+01
## countryUganda:occupationHomemaker 3.929e+01
## countryZimbabwe:occupationHomemaker 1.132e+01
## countryCameroon:mobileYes feature phone -1.014e+01
## countryIndia:mobileYes feature phone -1.203e+01
## countryIndonesia:mobileYes feature phone -6.720e+00
## countryKenya:mobileYes feature phone -1.023e+01
## countryNepal:mobileYes feature phone -8.874e+00
## countryPhilippines:mobileYes feature phone -1.090e+01
## countryUganda:mobileYes feature phone -1.256e+01
## countryZimbabwe:mobileYes feature phone -7.224e+00
## countryCameroon:mobileYes smart phone -7.591e+00
## countryIndia:mobileYes smart phone -8.314e+00
## countryIndonesia:mobileYes smart phone -6.244e+00
## countryKenya:mobileYes smart phone -8.372e+00
## countryNepal:mobileYes smart phone -6.467e+00
## countryPhilippines:mobileYes smart phone -1.035e+01
## countryUganda:mobileYes smart phone -9.429e+00
## countryZimbabwe:mobileYes smart phone -8.240e+00
## countryCameroon:internetYes -7.487e-02
## countryIndia:internetYes 1.909e-01
## countryIndonesia:internetYes 1.735e+01
## countryKenya:internetYes -1.775e-01
## countryNepal:internetYes 4.916e-01
## countryPhilippines:internetYes 1.163e+00
## countryUganda:internetYes 3.485e-01
## countryZimbabwe:internetYes -4.219e-01
## countryCameroon:educationPrimary school 2.102e+01
## countryIndia:educationPrimary school 8.260e+00
## countryIndonesia:educationPrimary school 2.197e+01
## countryKenya:educationPrimary school 1.578e-01
## countryNepal:educationPrimary school 6.556e+00
## countryPhilippines:educationPrimary school 4.262e+00
## countryUganda:educationPrimary school 1.154e+01
## countryZimbabwe:educationPrimary school 2.447e+01
## countryCameroon:educationSecondary school 1.169e+00
## countryIndia:educationSecondary school -9.693e+00
## countryIndonesia:educationSecondary school 1.699e+00
## countryKenya:educationSecondary school -9.395e+00
## countryNepal:educationSecondary school -9.995e+00
## countryPhilippines:educationSecondary school -9.008e+00
## countryUganda:educationSecondary school -1.141e+01
## countryZimbabwe:educationSecondary school 3.859e+00
## countryCameroon:educationDiploma 3.269e-01
## countryIndia:educationDiploma -1.120e+01
## countryIndonesia:educationDiploma 9.850e-01
## countryKenya:educationDiploma -1.107e+01
## countryNepal:educationDiploma -1.143e+01
## countryPhilippines:educationDiploma -1.079e+01
## countryUganda:educationDiploma -1.321e+01
## countryZimbabwe:educationDiploma 1.734e+00
## countryCameroon:educationDegree -2.979e-01
## countryIndia:educationDegree -1.286e+01
## countryIndonesia:educationDegree 2.862e-01
## countryKenya:educationDegree -1.025e+01
## countryNepal:educationDegree -1.340e+01
## countryPhilippines:educationDegree -1.190e+01
## countryUganda:educationDegree -9.794e+00
## countryZimbabwe:educationDegree 8.376e-01
##
## Intercepts:
## Value Std. Error t value
## 1|2 2.706380e+01 1.076900e+00 2.513160e+01
## 2|3 2.775100e+01 1.077600e+00 2.575330e+01
## 3|4 2.845470e+01 1.078200e+00 2.638990e+01
## 4|5 2.895350e+01 1.078700e+00 2.684100e+01
##
## Residual Deviance: 5866.142
## AIC: 6166.142
# compare with model 1
anova(m1, m2, test = "Chisq")
## Likelihood ratio tests of ordinal regression models
##
## Response: papm
## Model
## 1 country + gender + age + education + marital + occupation + live_area + mobile + internet
## 2 country + gender + age + education + marital + occupation + live_area + mobile + internet + gender:country + age:country + marital:country + occupation:country + mobile:country + internet:country + education:country
## Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 2864 6159.317
## 2 2744 5866.142 1 vs 2 120 293.1756 1.110223e-16
Anova(m2)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Analysis of Deviance Table (Type II tests)
##
## Response: papm
## LR Chisq Df Pr(>Chisq)
## country 1052.99 8 < 2.2e-16 ***
## gender 0.41 1 0.5230524
## age 3.39 1 0.0654129 .
## education 28.66 4 9.161e-06 ***
## marital 3.83 2 0.1470023
## occupation 7.84 4 0.0975951 .
## live_area 1.64 3 0.6506180
## mobile 15.12 2 0.0005211 ***
## internet 2.95 1 0.0860821 .
## country:gender 10.36 8 0.2405672
## country:age 22.71 8 0.0037536 **
## country:marital 31.38 16 0.0120169 *
## country:occupation 59.74 32 0.0020828 **
## country:mobile 41.60 16 0.0004521 ***
## country:internet 28.18 8 0.0004409 ***
## country:education 70.05 32 0.0001168 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# m2 is beter than m1.
https://peopleanalytics-regression-book.org/ord-reg.html:
dat_test1 <- dat1 %>% mutate(
test1 = ifelse(papm == "1", 0, 1),
test2 = ifelse(papm == "1" | papm == "2", 0, 1),
test3 = ifelse(papm == "1" | papm == "2" | papm == "3", 0, 1),
test4 = ifelse(papm == "1" | papm == "2" | papm == "3" | papm == "4", 0, 1)
)
table(dat_test1$papm)
##
## 1 2 3 4 5
## 1596 291 285 179 543
m1 <- glm(test1 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
data = dat_test1 %>% filter(country == "Bangladesh"),
family = "binomial"
)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
m2 <- glm(test2 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
data = dat_test1 %>% filter(country == "Bangladesh"),
family = "binomial"
)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
m3 <- glm(test3 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
data = dat_test1 %>% filter(country == "Bangladesh"),
family = "binomial"
)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
m4 <- glm(test4 ~ gender + age + education + marital + occupation + live_area + mobile + internet,
data = dat_test1 %>% filter(country == "Bangladesh"),
family = "binomial"
)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
coefficient_comparison <- data.frame(
test1 = summary(m1)$coefficients[, "Estimate"],
test2 = summary(m2)$coefficients[, "Estimate"],
test3 = summary(m3)$coefficients[, "Estimate"],
test4 = summary(m4)$coefficients[, "Estimate"],
diff1_2 = summary(m2)$coefficients[, "Estimate"] - summary(m1)$coefficients[, "Estimate"],
diff1_3 = summary(m3)$coefficients[, "Estimate"] - summary(m1)$coefficients[, "Estimate"],
diff1_4 = summary(m4)$coefficients[, "Estimate"] - summary(m1)$coefficients[, "Estimate"]
)
coefficient_comparison
## test1 test2 test3
## (Intercept) -38.43425101 -38.43425101 -40.93666714
## genderFemale 0.38700242 0.38700242 0.09338424
## age 0.08512152 0.08512152 0.08291879
## educationPrimary school 1.53180401 1.53180401 1.49493880
## educationSecondary school 17.77258036 17.77258036 18.95926562
## educationDiploma 19.23192113 19.23192113 19.17753807
## educationDegree 19.54783604 19.54783604 20.62539567
## maritalMarried/living with a partner -0.65031291 -0.65031291 -0.56711481
## maritalSeparated/widowed/divorced -2.11973287 -2.11973287 -1.11113159
## occupationUnemployed -18.58084035 -18.58084035 -19.58069143
## occupationStudent 0.11277406 0.11277406 -0.67586541
## occupationRetiree -20.44333003 -20.44333003 -20.16806042
## occupationHomemaker -16.21747376 -16.21747376 -16.88829039
## live_areaPeri-urban -0.45811434 -0.45811434 1.00093169
## live_areaRural 0.77632174 0.77632174 2.28880765
## mobileYes feature phone 14.16567476 14.16567476 14.22825219
## mobileYes smart phone 14.11045566 14.11045566 14.20706684
## internetYes 0.28136377 0.28136377 0.38110277
## test4 diff1_2 diff1_3
## (Intercept) -61.18243089 0 -2.502416133
## genderFemale -0.11672626 0 -0.293618176
## age 0.03502682 0 -0.002202736
## educationPrimary school 1.07945192 0 -0.036865204
## educationSecondary school 19.16571382 0 1.186685264
## educationDiploma -2.44021033 0 -0.054383058
## educationDegree 20.91176452 0 1.077559635
## maritalMarried/living with a partner 17.96932349 0 0.083198108
## maritalSeparated/widowed/divorced 18.74805329 0 1.008601275
## occupationUnemployed -16.15798504 0 -0.999851084
## occupationStudent -9.07947685 0 -0.788639469
## occupationRetiree -1.73019043 0 0.275269612
## occupationHomemaker -17.96936153 0 -0.670816628
## live_areaPeri-urban 18.80811464 0 1.459046033
## live_areaRural 20.42660607 0 1.512485918
## mobileYes feature phone -0.59054383 0 0.062577434
## mobileYes smart phone 5.33094271 0 0.096611175
## internetYes -5.93206007 0 0.099739006
## diff1_4
## (Intercept) -22.74817988
## genderFemale -0.50372868
## age -0.05009471
## educationPrimary school -0.45235209
## educationSecondary school 1.39313347
## educationDiploma -21.67213146
## educationDegree 1.36392848
## maritalMarried/living with a partner 18.61963640
## maritalSeparated/widowed/divorced 20.86778615
## occupationUnemployed 2.42285531
## occupationStudent -9.19225091
## occupationRetiree 18.71313961
## occupationHomemaker -1.75188777
## live_areaPeri-urban 19.26622898
## live_areaRural 19.65028434
## mobileYes feature phone -14.75621859
## mobileYes smart phone -8.77951296
## internetYes -6.21342384
We could see there is a large difference coefficients. The assumption of proportional odds is violated. We need to use different methods. The book recommend some other methods explain more details in Agresti (2010))
I still fit the model to individual country to see where the issues occur. When the absolute Value of the predictor is greater than 10, I make table for that independent variable and the outcome for that country. ## Calculate the confident interval I use Wald test to obtain the CI because the likelihood ratio test doesn’t work
countryList <- c("Bangladesh", "Cameroon", "India", "Indonesia", "Kenya", "Nepal", "Philippines", "Uganda", "Zimbabwe")
fun <- function(data, countryList) {
allData <- list()
for (ctry in countryList) {
m <- polr(formula = papm ~ gender + age + education + marital + occupation + live_area + mobile + internet, data = dat1 %>% filter(country == ctry), Hess = TRUE)
ci <- confint.default(m)
tmp <- as.data.frame(round(exp(cbind(OR = coef(m), ci)), 2)) %>%
mutate(
country = ctry,
predictor = row.names(m)
)
allData[[ctry]] <- tmp
}
combinedData <- do.call(rbind, allData)
return(combinedData)
}
tmp <- fun(dat1, countryList)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in polr(formula = papm ~ gender + age + education + marital + occupation
## + : design appears to be rank-deficient, so dropping some coefs
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in polr(formula = papm ~ gender + age + education + marital + occupation
## + : design appears to be rank-deficient, so dropping some coefs
## Warning in polr(formula = papm ~ gender + age + education + marital + occupation
## + : design appears to be rank-deficient, so dropping some coefs
tmp1 <- tmp %>%
mutate(independant_var = row.names(tmp)) %>%
separate(independant_var, c("count", "var"), sep = "\\.") %>%
mutate(var = str_replace_all(var, c(
"gender" = "gender.",
"age" = "age.Age",
"education" = "education.",
"marital" = "marital.",
"occupation" = "occupation.",
"live_area" = "live_area.",
"mobile" = "mobile.",
"internet" = "internet."
))) %>%
separate(var, c("var", "level"), "\\.") %>%
mutate(
OR = ifelse(OR >= 1000, 0, OR),
`2.5 %` = ifelse(`2.5 %` >= 1000, 0, `2.5 %`),
`97.5 %` = ifelse(`97.5 %` >= 1000, 0, `97.5 %`)
) %>%
mutate(level = as.factor(level))
tmp2 <- tmp1 %>%
filter(var == "gender")
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1")+
labs(title = "Reference - Male"
)
ggplotly(p)
tmp2 <- tmp1 %>%
filter(var == "age") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1")+
labs(title = "Reference - Contunious"
)
ggplotly(p)
tmp2 <- tmp1 %>%
filter(var == "education") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1")+
labs(title = "Reference - No formal education"
)
ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Warning: `position_dodge()` requires non-overlapping x intervals
tmp2 <- tmp1 %>%
filter(var == "marital") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1")+
labs(title = "Reference - Never married"
)
ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Warning: `position_dodge()` requires non-overlapping x intervals
tmp2 <- tmp1 %>%
filter(var == "occupation") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1") +
labs(title = "Reference - Employed"
)
ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Warning: `position_dodge()` requires non-overlapping x intervals
tmp2 <- tmp1 %>%
filter(var == "live_area") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1")+
labs(title = "Reference - Urban"
)
ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Warning: `position_dodge()` requires non-overlapping x intervals
tmp2 <- tmp1 %>%
filter(var == "mobile") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1") +
labs(title = "Reference - No"
)
ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Warning: `position_dodge()` requires non-overlapping x intervals
tmp2 <- tmp1 %>%
filter(var == "internet") %>%
mutate(level = as.factor(level))
p <- ggplot(data = tmp2, aes(x = level, y = OR, colour = country)) +
geom_line(group = 1, position = position_dodge(0.3)) +
geom_errorbar(aes(ymin = `2.5 %`, ymax = `97.5 %`), position = position_dodge(0.3), width = 2) +
geom_point(size = 3, position = position_dodge(0.3)) +
scale_y_continuous(trans = "log10") +
scale_fill_brewer(palette = "Set1") +
labs(title = "Reference - No"
)
ggplotly(p)
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis
We can use anova() function to obtain the P value for each independent variable using likeli-hood ratio test